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Riboswitches, structured elements in the untranslated regions of messenger RNAs, regulate gene 
expression by binding specific metabolities. We introduce a kinetic network model that describes 
the functions of riboswitches at the systems level. Using experimental data for flavin mononu- 
cleotide riboswitch as a guide we show that efficient function, implying a large dynamic range 
without compromising the requirement to suppress transcription, is determined by a balance be- 
tween the transcription speed, the folding and unfolding rates of the aptamer, and the binding rates 
of the metabolite. We also investigated the effect of negative feedback accounting for binding to 
metabolites, which are themselves the products of genes that are being regulated. For a range of 
transcription rates negative feedback suppresses gene expression by nearly 10 fold. Negative feed- 
back speeds the gene expression response time, and suppresses the change of steady state protein 
concentration by half relative to that without feedback, when there is a modest spike in DNA con- 
centration. A dynamic phase diagram expressed in terms of transcription speed, folding rates, and 
metabolite binding rates predicts different scenarios in riboswitch-mediated transcription regulation. 

PACS numbers: 



Introduction 



Riboswitches are ds-acting RNA elements located in 
the untranslated region of mRNAs that regulate associ- 
ated gene expression by sensing and binding target cellu- 
lar metabolites [T]-[3] . In bacteria, binding of metabolites 
to the conserved aptamer domain allosterically alters the 
folding patterns of the downstream expression platform, 
whose conformation controls transcription termination or 
translation initiation [3-5 . The target metabolites are 
usually the products or their derivatives of the down- 
stream gene that riboswitches control. Hence, metabo- 
lite binding to riboswitches serves as a feedback sig- 
nal to control RNA transcription or translation initia- 
tion. The feedback through metabolite binding is nat- 
urally designed to be a fundamental network motif for 
riboswitches. For example, tandem riboswitches respond 
to multiple metabolites to control a single gene with 
greater regulatory complexity [6j [7], while single glmS 
riboswitch has been shown to respond to multiple signals 
using both negative and positive feedback [8 j. Under- 
standing the various in vivo riboswitch functions requires 
a theoretical framework that takes into account the inter- 
play between speed of RNA transcription, folding kinetics 
of the nascent RNA transcript, and kinetics of metabo- 
lite binding to the nascent RNA transcript, and the role 
of feedback arising from interactions between synthesized 
metabolities and the transcript. The effects of speed of 
RNA transcription and metabolite binding kinetics have 
been examined in vitro in an insightful study involving 
the flavin mononucleotide (FMN) riboswitches [9]. They 
argued that FMN riboswitch is kinetically driven imply- 
ing that the riboswitch does not reach thermodynamic 
equilibrium with FMN before a decision between contin- 
ued transcription and transcription termination needs to 
be made. 



The regulatory roles played by riboswitches have also 
inspired design of novel RNA-based gene-control ele- 
ments that respond to small molecules [lOj [11] . Several 
models have been proposed to describe how riboswitches 
function and meet their regulatory demands p~2j [13] . 
However, they focused solely on the transcription pro- 
cess without accounting for the feedback effect from 
the metabolite produced by the gene encoding the ri- 
boswitch. Here, we introduce a general kinetic network 
model that can be used to describe both in vivo and in 
vitro functions of riboswitches. Our coarse-grained ki- 
netic network model, which takes into account the in- 
terplay of co-transcriptional folding, speed of transcrip- 
tion, and kinetics of metabolite binding, also models ef- 
fects of negative feedback loop so that predictions for in 
vivo functions of riboswitches can be made. As an illus- 
tration of the theory, we first consider the dependence 
of metabolite concentration on the regulation of in vitro 
transcription termination of FMN riboswitches without 
feedback loop, which enables us to obtain the range of 
folding rates and transcription rates that produce results 
consistent with experiments [9]. We then include the 
negative feedback loop in the network to study how ri- 
boswitches regulate gene expression at the systems level. 

General Kinetic Model 

The riboswitch is transcribed from the leader, the non- 
protein-coding region, of the associated gene (Fig. 1). 
We simplify the multiple complex in vitro biochemical 
steps in the function of the off riboswitch, involving tran- 
scription, co-transcriptional RNA folding and metabolite 
binding, to a few key kinetic steps (Fig. 1). With- 
out feedback, the first stage is the transcription of the 
aptamer domain (B). The antiterminator sequence is 
transcribed (B2) in the second step. At each stage, the 
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aptamer domain of the RNA transcript can either co- 
transcriptionally fold or unfold. Only when the aptamer 
domain is folded (£?*, can the RNA transcript bind 
the metabolite (M). At the second stage, when the ap- 
tamer domain is unfolded, the RNA transcript is in an 
alternative folding pattern with the formation of antiter- 
minator stem (B2). The final stage of the transcription 
occurs when the terminator sequence is transcribed (Ri). 
If the terminator sequence is transcribed following B 2 , 
the antiterminator structure prevents the formation of 
terminator stem and the transcription proceeds till the 
downstream coding region is fully transcribed (Rf). If 
the terminator sequence is transcribed following or 
£?|M, the absence of antiterminator allows the termina- 
tor to form, which subsequently leads to the dissociation 
of RNA transcript (B% t and B$ t M) from the DNA tem- 
plate and terminates the transcription. The feedback 
effect involved with translation and metabolite synthesis 
will be discussed in later sections. 

To assess how the metabolite concentration, [M], regu- 
lates transcription termination, we computed the fraction 
of terminated transcript, f ter , given an initial concen- 
tration of RNA transcript with aptamer sequence tran- 
scribed (B). Some of the rate constants can be estimated 
from the in vitro experiments [9] for FMN riboswitches, 
which we use to illustrate the efficacy of the theory. The 
experimental values of the FMN association rate constant 
kb for the FMN aptamer is ~ 0.1 /iM _1 s _1 , and the dis- 
sociation rate constant k-b is ~ 10 -3 s _1 , giving the 
equilibrium = k-b/^b = 10 nM. RNA polymerases 
(RNAP) pause at certain transcription sites during tran- 
scription. There are two pause sites for the FMN ri- 
boswitch, one after the aptamer domain sequence with a 
lifetime of the paused complex being about 10 s, and the 
other at the end of the antiterminator sequence with a 
lifetime ~ 1 min. To approximately account for the pause 
times in our simplified model, we observe that B2 repre- 
sents the transcript of the FMN riboswitch with part of 
the antiterminator out of RNAP, when RNAP pauses at 
the second pause site. Even with only part of the antiter- 
minator sequence, the transcript still has high probability 
of forming alternate folding patterns [9], similar to a full 
antiterminator sequence. Hence, we set the effective tran- 
scription rates k t \ = 0.1 s _1 and k t 2 = 0.016 s _1 , which 
reflects the pause times for the FMN ribsowtich (see Fig. 
IB for additional explanation for this approximation). 



riboswitch [9]. When the aptamer sequence is tran- 
scribed, the transcript favors the aptamer-folded state, 
and when the antiterminator sequence is transcribed, 
the folding pattern changes in favor of forming the an- 
titerminator stem with disruption of the aptamer folded 
structure [9]. Thus, there are restraints on the folding 
rates, kfi > and kf2 < &-/2- We also assume the 

same association (dissociation) rate constant for metabo- 
lite binding to B* , B\ , and B^ t because there is little 
change in the results when the values of kb2 and kbz are 
drastically altered. In addition, because we only allow 
the folded states of the aptamer to bind M the effective 
Kd is a convolution of the folding rates and the binding 
rate. Thus, even though kb is the same the decrease in 
the folding rate as the transcript length increases effec- 
tively decreases Kjj. The values of the transitions rates 
that reproduce the measured f ter (blue squares in Fig. 
2) are listed in Table I. The folding rates kfi and kfi 
are within an order of magnitude of the theoretical pre- 
diction based on, kf ~ k$e~^ \ ko ~ 10 6 s _1 , where 
N is the number of nucleotides [141 US] • Moreover, the 
rate kfi we obtained is in the same order of magnitude 
of the folding rate of other riboswitch aptamer with sim- 
ilar lengths observed in experiments [T6] . Thus, for the 
purposes of quantitatively describing the in vitro experi- 
ments we need only two kinetic rates (fc_/i and fc_/2)- 

The folding rate of the aptamer is comparable to k t i 
(Fig. 1). The transition rate (fc_/2) from B^ (Fig. 1) to 
B2 is 2-3 times the rate of transcription enlongation to 
the stage where terminator sequence is transcribed (£^2). 
Since the results are not sensitive to if the other 

rates are fixed, we choose fc_/i ~ &-/2 because both in- 
volve unfolding of the folded aptamer structure. With 
this assumption, the parameter set that emerge when 
our model is used to quantitatively describe (see Fig. 2) 
the in vitro kinetic experiments is unique. In addition, 
under these conditions, the regulation of transcription 
termination works when [M] is in large excess over RNA 
transcript, when [M] Q / [B] Q > 10. With [M] = 1 /iM 
for metabolite concentration and kb — 0.1 /iM _1 s _1 , the 
binding time is of ~10 s, which is of the same order of 
magnitude as the transcription elongation rate and the 
folding rate of the antiterminator stem. Consequently, 
the metabolite binding is unlikely to reach equilibrium 
before formation of the antiterminator stem. Large ex- 



Extraction of minimal set of parameters from in vitro 
transcription experiments 

To make testable predictions using our model, we need 
estimates of the co-transcriptional folding and unfold- 
ing rates of the aptamer B as well as £?2, the aptamer 
with the antiterminator sequence. The kinetic model de- 
scribed mathematically in the Supplementary Informa- 
tion (SI) can be used to extract parameters that most 
closely fit the measured dependence on [M] for the FMN 



kfi a k-fi kf 2 k-f 2 hi b k t 2 b kb b,c k-b b 



0.1 0.04 2.5xl0~ 3 0.04 0.1 0.016 0.1 10" 3 



a In unit of s _1 for all rates except k^. 
b Values from experimental data [9]- 
c In unit of uM~ 1 s~ 1 . 

TABLE I: Kinetic parameters for model in Figure 1 without 
feedback. 
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FIG. 1: (A) Kinetic network model for RNA transcription mediated by riboswitches. The leader, upstream of the protein- 
coding gene, consists of sequences that can be transcribed to the aptamer domain (£?), antiterminator (B2) and terminator 
region (Ri) of the riboswitch. After transcription initiation, elongation, folding of RNA transcript and metabolite binding 
are simplified to several key steps. Starting from the transcript B, where the aptamer sequence is transcribed, transcription 
can continue to B2 (antiterminator sequence transcribed) at a transcription rate constant kti. Further elongation through 
the terminator sequence with transcription rate constant kt2 results in the synthesis of full RNA without termination. Ri is 
the transcript withthe sequence of the protein-coding region starting to be transcribed, and eventually grows to Rf, the full 
protein-coding region transcribed, with a rate of k t 3- Besides transcription elongation, each of the transcript states, B and B 2 , 
can form states with aptamer domain folded (B* and B%) with a folding rate of kfi and respectively. The aptamer- folded 
states can bind metabolite (M) leading to the bound states (B*M and BJM) with association rate constant fe,. The transcripts 
in state B^ and B^ M can continue elongating till the terminator sequence is transcribed with their expression platform forming 
a transcription terminator stem and dissociate from the DNA template terminating transcription, with a rate of kter (B% t and 
B^M). The fraction of transcription termination, fter, is determined from the amount of the terminated transcripts (in green 
block) versus non-terminated transcripts (in blue block). In the presence of negative feedback loop (steps included by the 
box in dashes) additional biochemical steps have to be included. In this case after RNA is fully synthesized, it can produce 
protein P at a rate ki or get degraded with a rate kai- The fate of P is either degradation (rate kd2) or production of an 
inactive metabolite Mo, which is activated by the enzyme (E) encoded by the gene Of- The activated enzyme can bind to 
the folded aptamer and can abort transcription. (B) Simplification of the step from B2 to Ri for FMN riboswitch. In the 
application to FMN riboswitch, B2 represents the transcript out of the RNA polymerase at the second pause site [9]- The 
step B 2 — > Ri is a simplification of potential multiple chemical process, including pausing and emerging of the antiterminator 
sequence (B2 — > B3), and transcription to the terminator sequence (B3 — > Ri). The rate kt2 is approximated as the pausing 
rate k p , because pausing is likely to be the rate limiting step in the transcription process. 
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cess of metabolite, exceeding the equilibrium Ku (~ 10 
nM) for FMN binding to over RNA transcript is 
needed for the metabolites to bind the riboswitches with 
sufficient speed to regulate transcription under the con- 
ditions explored in experiments [9]. 



Results 
Dependence of fter on K2 

We investigated how f ter depends on variations in the 
transition rates around the parameter set listed in Ta- 
ble I. Figure 2 shows that the fraction of terminated 
transcripts converges to the same value in the limit of 
high metabolite concentration, independent of K 2 = 
k-f2/kf2 : while keeping the other paramters fixed. At 
high [M], B* and are always metabolite bound, which 
results in very low [B^]. Hence, varying fc_/2 does not 
affect fter at high [M] . In the limit of low [M] , f ter de- 
creases as K2 increases because fc_/2 exceeds the effective 
binding rate fc&[M] so that B2 is preferentially populated. 

The effective metabolite concentration T50, at which 
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where f drifter) * s the vame of fter at high (low) [M], 
does not change much as K2 is varied (see the inset in 
Fig. 2). Even when K2 is small T^q/Kjj > 1 implying 
the concentration of [M] is in excess of the equilibrium 
Kd to effect binding to B\. Since the population of B2 
is favored as K2 increases it follows that f ter decreases 
at all concentrations of the metabolite as K2 is changed 
from a low to a high value (Fig. 2). In addition, the tran- 
scription rate from B2 (or to the next stage where 
terminator sequence forms (fc^) is about one order of 
magnitude larger than kf2, which means that at low or 
normal metabolite concentration, transitions between B2 
and B2 states do not reach equilibrium before the termi- 
nator sequence is transcribed - a result that also follows 
from the inequality T 50 > K^. Finally, the dissociation 
rate constant is much smaller than k t 2-> which indicates 
that once the metabolite is bound, the bound state re- 
mains stable through transcription termination. Hence, 
the riboswitch is in the kinetically driven regime with the 
parameters used here - a conclusion that was reached in 
the previous study [9]. 



Dynamic Range and Thermodynamic Control 

Thermodynamic equilibrium between B2 and B^ can 
be reached only if the transcription speed is much slower 
than the transition rates between different folding pat- 
terns and the association rate with metabolites (Fig. 1). 
We varied the transcription speed to probe how the ri- 
boswitch can be driven from kinetic to thermodynamic 
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FIG. 2: Dependence of transcription termination on metabo- 
lite concentration without feedback. (A) Fraction of ter- 
minated RNA transcripts, fter, as a function of the loga- 
rithm of the metabolite concentration for different values of 
K2 = k-f2/kf2 with kb2 = kb3 = kb. Parameters that re- 
produce the in vitro experimental fter are listed in Table I. 
The inset in (A) shows the half-response metabolite concen- 
tration, T50, as a function of K2. (B) Sensitivity of fter to 
different values of kb2 and kb3- Except for modest changes in 
T50 there is little change in f ter when the binding rates are 
drastically altered. The blue points in (A) and (B) are data 
from experiments on FMN riboswitch [9]. 



control, which can be experimentally realized by increas- 
ing the pausing time, achievable by adding transcription 
factors, such as NusA. The dependence of f ter on [M] 
at various 72 = ^2/^/2 values shows that, in the limit 
of low metabolite concentration, f ter is roughly equal to 
the fraction of folded aptamer /#* ~ fcf 2/(^/2 + &-/2) ~ 
0.06. At high metabolite concentrations, almost all ri- 
boswitches are metabolite bound, and transcription is 
terminated with high probability (Fig. 3A). As 72 in- 
creases, the system transitons to a kinetically driven 
regime and the probability that the transcription is ter- 
minated at high [M] decreases, while the fraction of ter- 
minated transcript at low [M] increases. 

Interestingly, at high [M] we find that f ter decreases 
as 72 increases because in this limit the folded B^ has 
insufficient time to make a transition to B2. As a result 
the population of B2 decreases at high 72 resulting in a 
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FIG. 3: Speed of transcription and gene expression. (A) Fraction of terminated RNA transcripts, fter, as a function of the 
logarithm of metabolite concentration for different values of 72 = ^2/^/2- The parameters that reproduce the experimental 
fter results in 72 = 6. The inset shows log(Tso) as a function of log(72). (B) The dynamic range 77 = f^ r — fter, where f£> r 
{fter) is the value of f ter at high (low) metabolite concentration, as a function of K 2 . (C) Variation of T 50 as a function of K 2 . 
(D) Fraction of terminated RNA transcripts as a function of the logarithm of metabolite concentration for different values of 
K2 with 72 = 10 -2 . Other parameters used are listed in Table I. 



reduction in f ter (Fig. 3A) . Surprisingly, the exact oppo- 
site result is obtained at low [M] as 72 is varied. At low 
[M] and small 72 the binding rate fc& [M] is small enough 
that the transition to B2 occurs with high probability 
resulting in a decrease in f ter . As 72 increases the flux 
from B2 to B2 decreases, and the pathway to B\ from B* 
becomes relevant leading to an increase in f ter at low [M] 
(Fig. 3A). Thus, at high 72 and low [M] the extent of 
transcription termination is controlled by K\ = k-f\/kf\ 
and 71 = kn/kfi. The value of T50 increases substan- 
tially relative to Kjj as 72 increases (see the inset in Fig. 
3A). Even when 72 is very small T50 exceeds Kjj im- 
plying that it is difficult to drive the FMN riboswitch 
to thermodynamic control under the conditions used in 
experiments [9 . 

Figure 3B shows the dynamic range 77 = f^ r — fter as a 
function of K2 for different values of 72. The riboswitch 
functions with maximal dynamic range when the system 
is nearly under thermodynamic control corresponding to 
small K2 values. The range of T50 is between 0.1 \iM 
(w 100K D ) and 1 fiM (« 1000K D ) for 72 > 1 (Fig. 2C). 



When the unfolding rate of the aptamer folded struc- 
ture (#2) is comparable to the speed of transcription to 
the terminator sequence, K2 ~ 72, X50 has the small- 
est value. The minimum T50 decreases as 72 decreases, 
and T 50 becomes less dependent on K 2 when K 2 < 1. 
When 72 <C 1, as shown in Fig. 3D, the probability of 
transcription termination approaches unity in the limit 
of high metabolite concentration at all values of K2. On 
the other hand, in the limit of low metabolite concen- 
tration (fcb[M] is small), f ter increases and T50 decreases 
as K2 decreases. The results in Fig. 3 show that the 
efficiency of the riboswitch function is determined by a 
compromise between the need to maximize 77 (72 should 
be small) and the ability to terminate transcription (72 
should be large). 



Effect of aptamer folding rates on f ter 

In the kinetically driven regime (72 > 1), the probabil- 
ity of transcription termination depends on the fraction 
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FIG. 4: Aptamer folding rates and /* e r. (A) Fraction of 
terminated RNA transcripts as a function of log([M]) for dif- 
ferent values of K\ = k-fi/kfi using the parameters listed 
in Table I except for k-fi. (B) Fraction of terminated RNA 
transcripts as a function of log([M]) for different values of 
71 = kti/kfi. Inset shows the dependence of T50 on 71. 



of the aptamer-formed state before transcription of the 
antiterminator sequence. This fraction can be changed 
by altering fc/i, or equivalent ly the ratio K\ = k-fi/kfi, 
or by varying the transcription speed from the aptamer 
domain sequence to antiterminator sequence, k t \. When 
K\ ^> 1, most of the riboswitches do not form stable 
folded aptamer domain B* (Fig. 1), resulting in very 
small fraction of transcription termination and low re- 
sponse to changes in the metabolite concentration. As 
the transition rate from the unfolded state (B) to the 
aptamer-folded state (B*) increases relative to the re- 
verse transition rate, the fraction of terminated tran- 
scripts and dynamic range increases (Fig. 4A). In the 
limit of high [M] , the probability of transcription termi- 
nation approaches 1 when K\ <C 1, while in the low [M] 
limit, almost all the riboswitches are aptamer folded but 
without metabolite bound (B*) before transcription to 
antiterminator sequence. Just as in Fig. 2, T50 is not 
sensitive to changes in K\ (see inset in Fig. 4A). 

Figure 4B shows the concentration dependence of the 
fraction of terminated transcripts for different transcrip- 



tion rates to the antiterminator sequence, or the ratio 
7! = k t i/kfi. When the transcription rate is much 
faster than the folding rate, the riboswitch does not 
have enough time to form the aptamer-folded structure, 
which results in low fraction of terminated transcripts 
and low response to metabolite concentration change. 
When transcription rate is much slower than folding rate, 
the folded and unfolded states of the aptamer are able 
to reach equilibrium before transcription to antitermi- 
nator sequence. In the high metabolite concentration 
limit, the riboswitch is always metabolite bound result- 
ing in transcription termination. At low [M], the ri- 
boswitch does not bind the metabolite. The fraction 
of aptamer formed state without bound metabolite is 
kfi/(kfi + k-fi) ~ fter before transcription to the an- 
titerminator sequence. There is substantial variation in 
T 50 as 71 changes as the inset in Fig. 4B shows. Thus, 
besides the speed of transcription and the binding rates, 
the folding rates of the aptamers have considerable influ- 
ence on fter (compare insets in Fig. 3 A and Fig. 4B). 



Transcription with negative feedback loop 

Most riboswitches regulate gene expression of the 
downstream platform that encodes for proteins involved 
in the production of the specific metabolite that itself 
binds to the riboswitch (Fig. 1). Therefore, sensing and 
binding of its own metabolite by the riboswitch acts as 
a feedback to control gene expression. For riboswitches 
that suppress gene expression by binding to metabolites 
with high selectivity (for example, guanine riboswitches 
or FMN riboswitches) such feedback loop is an example 
of negative autoregulation, which has been widely stud- 
ied in gene regulation networks associated with transcrip- 
tion factors [2T] , We include the role negative feedback 
plays in regulating transcription termination by general- 
izing the in vitro kinetic model considered in the previous 
section. Our minimal model, illustrated in Fig. 1 and 
described in detail in the SI, provides a framework for 
interpreting future in vivo experiments. 

We consider transcription and translation in a cell and 
take into account RNA degradation and cell expansion. 
The transcription process is similar to that described 
without feedback loop, except now we include the effect 
of cell expansion and RNA degradation. We assume that 
the cell grows at a rate of \i = 5 x 10 -4 s _1 , resulting 
in a typical doubling time of ln2//i ~ 20 minutes for an 
E. coli cell, and that the degradation rate of the fully 
transcribed RNA or terminated RNA transcript is kd\. 
The values of fc^i? an d other parameters in the feedback 
loop are in Table II. The fully transcribed RNA serves 
as a template for the translation of the protein (P) that 
synthesizes metabolite Mo, which is then converted to an 
active form M by the enzyme E encoded by the gene Of 
(Fig. 1). The species M, with a degradation rate of ^3, 
is the target metabolite that binds to the riboswitch. 

In the case of FMN riboswitches, Mq represents ri- 
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d Ref.[T9l 

e Ref.[20] 

^See SI." 



TABLE II: Additional kinetic parameters for model in Figure 1 with negative feedback. 



boflavin, the eventual product encoded by the gene ribD 
(Or), that is subsequently converted to FMN by flavoki- 
nase (E), synthesized by the gene ribC (Of)- The degra- 
dation rate kd3 takes into account the effect of conversion 
from FMN to FAD (flavin adenine dinucleotide) by FAD 
synthetase in vivo (see SI for more details). However, we 
neglect the potential binding of FAD to the riboswitch 
because there is a 60 fold difference (or potentially even 
larger factor in the absence of FMN) in the binding of 
FMN and FAD to the FMN riboswitch [22 . In the model 
with negative feedback the extent of regulation by ri- 
boswitches is expressed in terms of the production of the 
protein P (Fig. 1). 

We assume that the activated level of the operon, Or, 
for transcription initiation of the riboswitch is a constant, 
and set [DN A] = [Or] = 2.5 nM, which is equivalent 
to one DNA molecule in an E. coli cell, and assume 
that the aptamer sequence, B, is produced with an ef- 
fective rate constant k\ — 0.016 s 1 , taking into account 
the transcription initiation rate (~ 50 s between initia- 
tion events; [23l[24]) and the typical transcription speed, 
~ 10 — 35 nucleotides per second [25] . With these param- 
eters fixed, which are used for illustrative purposes only, 
we can study the effect of feedback by varying the effec- 
tive rate kp, at which E is produced from Op. We set 
[Op] = [Or]. If E is produced at a rate similar to that 
of protein P without feedback, then kp ~ 1 s _1 , which 
we set as a reference rate kpo (see SI for details). At this 
rate, the steady state level of enzyme E is ~ 10 3 copies 
per cell. The variation of rate kp can result from delays 
or speed up in the process of the transcription from Of 
or translation of E, or even deficiency in Op. 

Dependence of [P] on enzyme production and 
metabolite binding rates 
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FIG. 5: Effect of negative feedback. (A) Protein (P in Fig. 
1) concentration at different 72 values as a function of the 
logarithm of fop, the production rate constant of enzyme E 
that produces the metabolite M, relative to /cfo = 1 s _1 . The 
parameters are given in Table I and II except k t 2- (B) The 
extent of regulation expressed in terms of protein level as a 
function of the logarithm of association rate constant kb for 
metabolite binding. The parameter are listed in Table I and 
II except kb and k-b- Kd is fixed at 10 nM. 



We assess the extent of regulation due to feedback by 
the changes in the protein level expression, [P], as the 
parameters in the network are varied. The results in 
Fig. 5 A show that when kp/kpo is low, very few active 
metabolites are formed to suppress protein expression. 
Consequently, the expression level of protein does not 
change at low kp/kpo at all values of 72 (Fig. 5A). This 
finding explains the observation that deficiency in ribC 



(Of in Fig. 1), the gene that encodes for flavokinase 
(E in Fig. 1), causes accumulation of riboflavin (Mq in 
Fig. 1) without converting it to FMN, and thus cannot 
suppress the synthesis of riboflavin [26j[27]. When kp in- 
creases to about 10 _3 /cfo, the expression of the protein 
starts to be suppressed. There is a substantial suppres- 
sion of protein concentration (Fig. 5A) at all 72 values 
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when kp/kpo > 10 -2 . As kp exceeds O.lkpo, the sup- 
pression begins to saturate (Fig. 5A) because most of 
Mq produced are converted to M, and the [P] level is 
essentially determined by the production of [M]o, which 
is independent of kp. Thus, at all 72 values, [P] level 
varies between two steady state values as kp is changed. 

In contrast to the results in Fig. 5A, the extent of 
regulation ([P] levels) varies greatly with binding rate 
constant of metabolites while keeping Kp> fixed at 10 
nM (Fig. 5B). Changes in [P] as fc& is varied, which af- 
fect synthesis of M (Fig. 1), depend on kp/kpo. When 
kp > kpo, most of the aptamer folded structures are 
metabolite bound at experimental value of fc& = 0.1 
/iM _1 s _1 , and hence the level of protein expression is 
suppressed. Therefore, with ~ 10 3 copies of enzyme E in 
a cell, the production of [P] decreases substantially even 
if the binding rate constant is small. Binding occurs be- 
cause the concentration of active metabolites (~ 25 jaM) 
is far larger than RNA transcripts (~ 10 nM), result- 
ing in a high effective binding rate fc&[M]. The ability 
to suppress protein production decreases if the binding 
rate constant is smaller than the experimental value by 
more than one order of magnitude, or when the value of 
kp decreases. Not surprisingly, when kp is very low, the 
dependence of binding rate on expression of P decreases. 
As a consequence the changes in P production decreases 
greatly as kp/kpo decreases (Fig. 5B). Thus, only over a 
small range of fc& and kp/kpo does the riboswitch func- 
tion with sufficient dynamic range. 



Interplay between co-transcriptional folding and 
transcription speed 
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FIG. 6: Role of kti on negative feedback. (A) The frac- 
tion of fully transcribed RNA, f tra = 1 — fter, as a func- 
tion of ktilk-ji. The numbers give values of /cp/fcpo with 
kFo = 1 s _1 . The dependence of ftra on ktijk-ji shows 
three regimes: thermodynamic controlled regime for low kti 
(kt2/k-f2 < 0.025, or kti < 0.1k-b), aptamer folding dom- 
inated regime for high kti (kti > k-fi), and intermediate 
regime with significant metabolite binding. (B) The P con- 
centration as a function of k t i/k-fi at different /cf/&fo values 
(as described in (A)). 



The dependence of transcription rate on the extent of 
regulation, shown in Fig. 6, exhibits three distinct func- 
tional modality depending on the value of k t 2 relative to 
k-f2- First, consider the case with kp = (black line in 
Fig. 6A). If k t 2 ^> k-f2, the fraction of fully transcribed 
RNA (Fig. 6A), f tra = 1 - f ter , becomes 



[RNA] + 



[RNA] 



V + k -n ) 



(2) 



where [RNA] is the sum of concentration of the fully 
transcribed and terminated RNA transcripts, and K\ = 
k-fi/kfi. In this limit, f tra depends predominantly on 
the folding transition rate before the antiterminator se- 
quence (B 2 ) is transcribed (Fig. 1). Hence, f tra is a func- 
tion of kti, and the rate of cell expansion. When 
k t 2/k-f2 decreases, co-transcriptional folding results in 
the formation of the antiterminator stem, which prevents 
transcription termination. If fc t2 <^ &-/2 an d [i <C &-/2 5 
then 



K 2 



/2 



K 2 + 1 k f2 



-/2 



(3) 



where K 2 = fc_/ 2 /%2- In this limit, B 2 and B\ (Fig. 1) 
are in equilibrium. These results are the same as those 
in the limit of low metabolite concentration without the 
feedback loop (Fig. 3A). 

For finite kp, when k t 2 is fast relative to fc_/2 the 
expression level of protein P is nearly independent of 
k t 2 (Fig- 6B). In this regime transcription termination 
and hence the extent of completed transcription is deter- 
mined by the folding rates of the aptamer, which is not 
greatly affected by negative feedback. When the tran- 
scription speed decreases, the expression of protein P in- 
creases for small kp (Fig. 6B). The expression level of P 
reaches a peak when ^2/^-/2 ~ 0.1 — 1, and it disappears 
when kp/kpo > 1 because the metabolite binding be- 
comes significant enough to stabilize the folded aptamer 
structure and offset the effect of formation of the an- 
titerminator stem. This is illustrated in Fig. 6A, which 
shows that f tra decreases significantly when k t 2 goes be- 
low 0.1/c_/2 and reaches a minimum at ^2/^-/2 ~ 0.01. 
The dependence of [P] on ^2/^-/2 is maximal when 
k t 2/k-f2 ~ 0.01 — 0.1. When the transcription rate de- 
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creases further (fc t2 < 0.1fc_b, or ^2/^-/2 < 2.5 x 10 -3 
as shown by the left dashed line in Fig. 6A), the fraction 
of fully transcribed RNA increases sharply because the 
transcription rate is slow enough for the dissociation of 
metabolites from riboswitches to occur significantly. The 
system can establish thermodynamic equilibrium, which 
increases the favorability of the formation of antitermi- 
nator stem when k t 2 decreases resulting in decreasing 
effective binding rate k^M]. However, the overall ex- 
pression level of protein P becomes very low because 
slow transcription results in a decrease in P production. 
In addition, there is also significant probability of RNA 
degradation, which also results in a decrease in P ex- 
pression. Therefore, the extent of regulation due to the 
negative feedback of metabolite binding has a maximal 
effect when ^2/^-/2 ~ 0.01 where the protein expression 
is suppressed by metabolite binding by as large as 85%. 



Dynamic phase diagram: Competition between 
folding, transcription, and binding rates 

To have a more complete picture on how the inter- 
play between folding of RNA transcripts, transcription, 
and metabolite binding regulate the expression of P, we 
study the dependence of [P] on the transcription rates 
and the effective binding rate k b [M] . The dynamic phase 
diagram in Fig. 7 is calculated by varying both k t \ (k t 2) 
and with Kb = 10 nM. In the first stage of transcrip- 
tion elongation, i.e. after the aptamer sequence is tran- 
scribed, the formation of the aptamer structure is the key 
step in regulating transcription termination. Thus, the 
folding rate kji and the effective metabolite binding rate 
are the key rates in competition with k t \ for regulation 
of [P]. Figure 7A shows three regimes for the depen- 
dence of [P] on k t i and k^M). In regime I, k t \ > fc/i, 
the folding rate is slow relative to transcription to the 
next stage (Fig. 1). Thus, the aptamer structure does 
not have enough time to form. The dominant flux is 
from B to £?2, which leads to high probability of fully 
transcribed RNA downstream because of the low tran- 
sition rate from B2 to B%. The metabolite binding has 
little effect on protein expression in this regime, partic- 
ularly for large fcti/fc/i, and hence the protein is highly 
expressed. In regime II, k^M] < k t \ < fc/i, the aptamer 
has enough time to fold but the metabolite binding is 
slow. The dominant flux is B — » B* — >> B%, leading to 
formation of antitermination stem (B% —> B2) or tran- 
scription termination — » B^). The expression level 
of protein is thus mainly determined by fc_/2 and fc^? 
and the protein production is partially suppressed in this 
regime. In regime III, k t i < fc/i and k t i < fc&[M], the 
aptamer has sufficient time to both fold and bind metabo- 
lite, the dominant pathway is B ->> B* ->> B*M ->> B\M, 
leading to transcription termination. The protein pro- 
duction is highly suppressed in this regime. The results 
using parameters from Table I and Table II (k t \ = kfi 
and kb[M]/kfi ~ 25) fall on the interface of regime I 



and regime III, as shown by the arrow in Fig. 7A. The 
metabolite binding does not reach thermodynamic equi- 
librium due to low dissociation constant. However, the 
effective binding rate is high because the steady state 
concentration of metabolites (~ 25 /iM) is in large ex- 
cess over RNA transcripts. Thus, the riboswitch is ki- 
netically driven under this condition even when feedback 
is included. 

With kti comparable to kfi, at the second stage of 
transcription elongation the key step against transcrip- 
tion termination is the formation of the antiterminator 
stem [B2 — > B2). Figure 7B also shows three regimes 
for the dependence of [P] on k t 2 and fc&[M] relative to 
fc-/2, and the associated most probable pathways are 
displayed on the right. In regime I, kb[M] > fc-/2, the 
effective binding rate competes favorably with — >> B2 
transition so that £?2, if populated, is not likely to form 
the antiterminator stem. However, in this regime, the 
effective binding rate is also likely to be larger than k t i, 
resulting in most of the metabolite binding occurring at 
the first stage of transcription. Protein production is 
partially suppressed with the flux towards transcription 
termination flowing through B*M — >> B%M — >> T (ter- 
minated transcript). In regime II, fc&[M] < k-f2 < k t 2, 
both the metabolite binding and — » B2 transition are 
too slow to occur. The protein production level is mainly 
determined by k t \ and kf\. The major pathways are 
B -> B 2 -> RNA and B -> B% -> T leading to 

partial protein suppression. In regime III, kb[M] < fc_/2 
and k t 2 < &-/2, metabolite binding is too slow to oc- 
cur, but the riboswitch has enough time to form the an- 
titerminator stem before the terminator sequence is tran- 
scribed. The major flux from B\ flows to B 2l leading to 
fully transcribed RNA, and proteins are highly expressed. 
We note that using the parameters from Table I and Ta- 
ble II (k t2 /k-f2 = 0.4 and k b M/k- f2 ~ 60), the result 
falls in regime I with partial protein suppresion. Among 
the three regimes, regime I has efficient negative feed- 
back, while the slow metabolite binding in regime II and 
regime III makes results resemble to those without feed- 
back. The dynamic phase diagrams predicts results with 
limiting cases of various parameters, whose values may 
be within range in vivo and most certainly in vitro. 



Role of feedback in response to DNA bursts 

To assess how feedback affects the response to a sud- 
den burst in DNA concentration we calculated the time 
dependent changes in the protein concentration, 

AP(t) = P[DNA] f (t) ~ P[DNA]M (4) 
^ r st [DNA]f [DNA]i 

when the DNA concentration is switched from [DNA]i 
to [DNA)f. In Eq. (4) st stands for steady state. Fig 
8 A shows the response time, defined as the time needed 
to reach half-way to the new steady-state level (dashed 
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FIG. 7: Dependence of protein production on the network parameters with feedback. (A) Protein levels as functions of 
hi/kfi and kb with negative feedback using parameters in Tables I and II. The scale for the [P] production is shown in the 
color spectrum. The dependence of [P] on k t i and kb is catagorized into three regimes (see text for details). Points on the 
dashed line separating regime II and regime III satisfy kb[M] — kti. The major pathway in thetranscription process in each 
regime is shown on the right. The arrow indicates the data point resulting from using the value of kti and kb in Table I. (B) 
Expression level of proteins as functions of kti/k-fi and kb with negative feedback using the parameters in Tables I and II. 
The dependence of [P] on kti and kb[M] is catagorized into three regimes. Points on the dashed line separating regime I and 
regime II /III satisfy kb[M] = k-fi. The corresponding major transcription pathways are shown on the right. The data point 
corresponding to the arrow results from using the value of kti and kb in Table I. 



line in Fig. 8A), for [DNA] f = 3.75 nM and [DNA]i 
= 2.5 nM using folding and unfolding rates one order of 
magnitude larger than those in Table I. Small transient 
fluctuations in DNA concentration could arise from envi- 
ronmental stresses, and hence it is interesting to examine 
the response of the network to such changes. The values 
of the folding rates from Table I results in little differ- 
ence in response time between cases with negative feed- 
back and without feedback. However, with larger folding 
rates, the response time for the systems with negative 
feedback is significantly shorter than without feedback 
(Fig. 8A). The fractional change of the fully transcribed 
RNA, AF(t) = f tra ([DNA] f ,t) - f tr a([DNA]i,t) (Fig. 
8B), shows a slight increase with overshoot initially be- 



fore settling into a steady state level lower than the orig- 
inal one in the case of negative feedback (blue line in Fig. 
8B). For the case without feedback, the fraction of fully 
transcribed RNA decreases initially before reaching the 
expected steady state level (red line in Fig. 8B). 

With negative feedback and 10- fold increase in over- 
all folding and unfolding rates, the fractional increase 
in the protein steady state level, 5P st = (P^dna] ~ 
P[DNA]i)/P[DNA]ii m res P onse to the increase in DNA 
level, SD = ([DNA]f - [DNA]i)/[DNA] u is reduced by 
more than half of that in the case without feedback for 
certain range of parameters, as shown in Fig. 8C and 
8D. Without feedback, SP st = SD. Negative feedback 



11 




log(k t2 /k_ f2 ) log(k t2 /k_ f2 ) 

FIG. 8: Role of feedback to a spike in the DNA concentration. (A) Response of protein level (AP) relative to the change of 
steady state level (AP st ) when the DNA level increases by 50%. The response time is 850 s with negative feedback (blue line) 
and 1450 s without feedback (red line). The values of folding and unfolding rates used are 10 fold larger than those given in 
Table I. (B) Response of the fraction of fully transcribed RNA, AF, when the DNA level increases by 50%. The fraction of 
fully transcribed RNA increases initially with negative feedback (blue line) and decreases initially without feedback (red line) 
before settling to new steady state. (C) Fractional change in protein steady state level relative to that in DNA concentration, 
5P st /SD, in repsonse to a 50% rise in DNA level, as a function of £^2/^-/2 and kb for the riboswitch network with negative 
feedback using paramters in Table I. Points on the dashed line separating regime I and regime II /III satifsy kb[M] = fc-/2- (D) 
Same as (C) with the overall folding and unfolding rates being 10 fold greater than those in Table I. The arrow indicates the 
data point when k t 2 and kb are at the values from Table I. 



noticeably reduces the variations of expression in protein 
due to DNA level change. Substantial reduction occurs 
when the effective binding rate is comparable to k_f2 
and when k t 2 < &-/2(the interface between regime I and 
regime III in Fig. 8D). 

Discussion 

Transcription, regulated by metabolite binding to ri- 
boswitches, depends on an interplay of a number of time 
scales that are intrinsic to the co-transcriptional folding 
of the riboswitch as well as those determined by cellullar 
conditions. For a riboswitch to function with a large 
dynamic range, transcription levels should change signif- 
icantly as the metabolite concentration increases from a 
low to high value. In the high concentration limit, RNA 
transcript in the aptamer folded state binds a metabolite. 



Low dissociation rate constant results in the formation 
of terminator stem, which subsequently terminates tran- 
scription. In the low concentration limit, the aptamer 
folded state is mostly unbound and can remain folded till 
transcription termination or can fold to the antitermina- 
tor state leading to the transcription of the full RNA. 
The levels of transcription termination is thus controlled 
by the transition rates between the aptamer folded and 
unfolded states. 

For in vitro description the efficiency of riboswitch is 
determined by two conflicting requirements. If 77, the dy- 
namic range, is to be maximized, then 72 has to be suf- 
ficiently low. However, at low 72 and realistic values of 
the metabolite concentration, f ter « 1 which implies the 
switching function (needed to abort transcription) can- 
not be achieved. Thus, 72 has to have an optimal range 
(72 ~ (1 — 10)) for the riboswitch to have sufficient dy- 
namic range without compromising the ability to switch 
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from an on to off state. 

In the presence of negative feedback loop the con- 
centration of target metabolites is also regulated by 
gene expression. Under nominal operating conditions 
(k t 2/k-f2 ~ 0.01 — 0.1) binding of target metabolites, 
products of the downstream gene that riboswitches reg- 
ulate, significantly suppresses the expression of proteins. 
Negative feedback suppresses the protein level by about 
half relative to the case without feedback. In vivo, the 
presence of RNA binding proteins, such as NusA, may 
increase the pausing times, thus effectively reducing the 
transcription rates. Thus, the repression of protein level 
by the riboswitch through metabolite binding may be up 
to 10 fold. Faster RNA folding and unfolding rates than 
those we obtained may also increase the suppression by 
negative feedback and broaden the range of transcription 
rates over which maximal suppression occurs. These pre- 
dictions are amenable to experimental test. 

In response to changes in the active operon level, the 
negative feedback speeds up the response time of ex- 
pression and modestly reduces the percentage change in 
the protein level relative to change in operon level. The 
steady state level of expression for autoregulation varies 
as a square-root of the DNA concentration. Adaptive bi- 
ological systems may minimize the variation in gene ex- 



pression to keep the systems functioning normally even 
when the environments change drastically. One may need 
to consider more complex networks than the single au- 
toregulation in the transcription network to find near per- 
fect adaptation to the environmental change [28] . 

Riboswitches provide novel ways to engineer biologi- 
cal circuits to control gene expression by binding small 
molecules. As found in tandem riboswitches [6j [7], multi- 
ple riboswitches can be engineered to control a single gene 
with greater regulatory complexity or increase in the dy- 
namic range of gene control. Synthetic riboswitches have 
been successfully used to control the chemotaxis of bac- 
teria [29]. Our study provides a physical basis for not 
only analyzing future experiments but also in anticipat- 
ing their outcomes. 
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